Ambulatory or residential? a multi-state analysis of treatments for substance use disorders (Step 1.5)

Parametric

Assessment of Fit of Transitions

We need to derive appropriate functional forms and define respective survival functions. One reason to favor patametric models is that they can be easier to generalize. Several candidate distributions were considered to model transitions, including Weibull (assume a monotonically increasing or decreasing hazard), Log-logistic (non-monotonic hazards), Gompertz (assume a monotonically increasing or decreasing hazard), Log-normal (non-monotonic hazards), Exponential (constant hazard), Gamma, Generalized gamma & Generalized F (more flexible).

The following plots fitted survival curves from each model (colored lines), with the Kaplan-Meier estimate, in red, obtained from an example of Jackson available here, added to the contributions of Wathers & Cutler available here.

#options(warn=-1)

n_iter<-10000

tiempo_antes_fits<-Sys.time()

#Weathers, Brandon and Cutler, Richard Dr., "Comparison of Survival Curves Between Cox Proportional 
#Hazards, Random Forests, and Conditional Inference Forests in Survival Analysis" (2017). All Graduate 
#Plan B and other Reports. 927. 
#https://digitalcommons.usu.edu/gradreports/927 

#<div style="border: 1px solid #ddd; padding: 5px; overflow-y: scroll; height:350px; overflow-x: scroll; width:100%">            
#https://devinincerti.com/2019/01/01/sim-mstate.html

n_trans <- max(trans_matrix, na.rm = TRUE)
fits_wei <- vector(mode = "list", length = n_trans)
fits_weiph <- vector(mode = "list", length = n_trans)
fits_llogis <- vector(mode = "list", length = n_trans)
fits_gomp <- vector(mode = "list", length = n_trans)
fits_logn <- vector(mode = "list", length = n_trans)
fits_exp <- vector(mode = "list", length = n_trans)
fits_gam <- vector(mode = "list", length = n_trans)
fits_ggam <- vector(mode = "list", length = n_trans)
fits_genf <- vector(mode = "list", length = n_trans)
fits_genf_orig <- vector(mode = "list", length = n_trans)
fits_ggam_orig <- vector(mode = "list", length = n_trans)
fits_rp02 <- vector(mode = "list", length = n_trans)
fits_rp03 <- vector(mode = "list", length = n_trans)
fits_rp04 <- vector(mode = "list", length = n_trans)
fits_rp05 <- vector(mode = "list", length = n_trans)
fits_rp06 <- vector(mode = "list", length = n_trans)
fits_rp07 <- vector(mode = "list", length = n_trans)
fits_rp08 <- vector(mode = "list", length = n_trans)
fits_rp09 <- vector(mode = "list", length = n_trans)
fits_rp10 <- vector(mode = "list", length = n_trans)
fits_c_wei <- vector(mode = "list", length = n_trans)
fits_c_weiph <- vector(mode = "list", length = n_trans)
fits_c_llogis <- vector(mode = "list", length = n_trans)
fits_c_gomp <- vector(mode = "list", length = n_trans)
fits_c_logn <- vector(mode = "list", length = n_trans)
fits_c_exp <- vector(mode = "list", length = n_trans)
fits_c_gam <- vector(mode = "list", length = n_trans)
fits_c_ggam <- vector(mode = "list", length = n_trans)
fits_c_genf <- vector(mode = "list", length = n_trans)
fits_c_genf_orig <- vector(mode = "list", length = n_trans)
fits_c_ggam_orig <- vector(mode = "list", length = n_trans)
fits_c_rp02 <- vector(mode = "list", length = n_trans)
fits_c_rp03 <- vector(mode = "list", length = n_trans)
fits_c_rp04 <- vector(mode = "list", length = n_trans)
fits_c_rp05 <- vector(mode = "list", length = n_trans)
fits_c_rp06 <- vector(mode = "list", length = n_trans)
fits_c_rp07 <- vector(mode = "list", length = n_trans)
fits_c_rp08 <- vector(mode = "list", length = n_trans)
fits_c_rp09 <- vector(mode = "list", length = n_trans)
fits_c_rp10 <- vector(mode = "list", length = n_trans)
km.lc <-vector(mode = "list", length = n_trans)

fits_cox<-list()
fits_c_cox<-list()
#"gengamma" Generalized gamma (stable parameterisation)
#"gengamma.orig"    Generalized gamma (original parameterisation)
#"genf" Generalized F (stable parameterisation)
#"genf.orig"    Generalized F (original parameterisation)
#"weibull"  Weibull
#"gamma"    Gamma
#"exp"  Exponential
#"lnorm"    Log-normal
#"gompertz" Gompertz

library(flexsurv)

#Specify the formula
fitform <- Surv(time, status) ~ 1

for (i in 1:n_trans){  
  fits_wei[[i]] <- flexsurvreg(formula=fitform,
                               data = subset(ms_d_match_surv, trans == i),
                               dist = "weibull")
}

for (i in 1:n_trans){  
  fits_weiph[[i]] <- flexsurvreg(formula=fitform,
                               data = subset(ms_d_match_surv, trans == i),
                               dist = "weibullph")
}

for (i in 1:n_trans){
  fits_llogis[[i]] <- flexsurvreg(formula=fitform,
                               data = subset(ms_d_match_surv, trans == i),
                               dist = "llogis")
}

for (i in 1:n_trans){
  fits_gam[[i]] <- flexsurvreg(formula=fitform,
                               data = subset(ms_d_match_surv, trans == i),
                               dist = "gamma")
}
#LOS ERRORES OCURREN CUANDO NO HAY COVARIABLES
#In (function (q, shape, rate = 1, scale = 1/rate, lower.tail = TRUE,  ... :   NaNs produced
#gamma no funcionó: NaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs produced
for (i in 1:n_trans){
  fits_ggam[[i]] <- flexsurvreg(formula=fitform,
                               data = subset(ms_d_match_surv, trans == i),
                               dist = "gengamma")
}

for (i in 1:n_trans){
  fits_gomp[[i]] <- flexsurvreg(formula=fitform,
                                data = subset(ms_d_match_surv, trans == i),
                                dist = "gompertz")
}

## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193: valor inicial en 'vmmin' no es finito

for (i in 1:n_trans){
  fits_logn[[i]] <- flexsurvreg(formula=fitform,
                                data = subset(ms_d_match_surv, trans == i),
                                dist = "lnorm")
}

for (i in 1:n_trans){
  fits_exp[[i]] <- flexsurvreg(formula=fitform,
                               data = subset(ms_d_match_surv, trans == i),
                               dist = "exp")
}

no_attempts <- 20

for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
          r <- flexsurvreg(formula=fitform,
                               data = subset(ms_d_match_surv, trans == i),
                               dist = "gengamma.orig")
      )
    }
    fits_ggam_orig[[i]] <- r
}  

for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
          r <- flexsurvreg(formula=fitform,
                               data = subset(ms_d_match_surv, trans == i),
                               dist = "genf.orig")
      )
    }
    fits_genf_orig[[i]] <- r
}  

## Warning in flexsurvreg(formula = fitform, data = subset(ms_d_match_surv, :
## Optimisation has probably not converged to the maximum likelihood - Hessian is
## not positive definite.

## Warning in flexsurvreg(formula = fitform, data = subset(ms_d_match_surv, :
## Optimisation has probably not converged to the maximum likelihood - Hessian is
## not positive definite.

for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <-  flexsurvreg(formula=fitform,
                               data = subset(ms_d_match_surv, trans == i),
                               dist = "genf")
      )
    }
    fits_genf[[i]] <- r
}  

for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform,k=2,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_rp02[[i]] <- r
}  
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform,k=3,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_rp03[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform,k=4,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_rp04[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform,k=5,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_rp05[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform,k=6,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_rp06[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform,k=7,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_rp07[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform,k=8,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_rp08[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform,k=9,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_rp09[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform,k=10,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_rp10[[i]] <- r
}

for (i in 1:n_trans){
km.lc[[i]] <- survfit(formula= fitform, data = subset(ms_d_match_surv, trans == i))
}

transition_label<- plots$title
attr(transition_label,"names")<- plots$trans


get_distinct_hues <- function(ncolor,s=0.5,v=0.95,seed=2125) {
  golden_ratio_conjugate <- 0.618033988749895
  set.seed(seed)
  h <- runif(1)
  H <- vector("numeric",ncolor)
  for(i in seq_len(ncolor)) {
    h <- (h + golden_ratio_conjugate) %% 1
    H[i] <- h
  }
  hsv(H,s=s,v=v)
}

distinct_hues<-get_distinct_hues(20)

layout(matrix(1:n_trans, nc = 2, byrow = F))
for (i in 1:n_trans){
plot(km.lc[[i]], col="red", conf.int=F);
  lines(fits_wei[[i]], col=distinct_hues[1], ci=F);
  lines(fits_weiph[[i]], col=distinct_hues[2], ci=F);
  lines(fits_gomp[[i]], col=distinct_hues[3], ci=F);
  lines(fits_llogis[[i]], col=distinct_hues[4], ci=F);#A0A36D
  lines(fits_gam[[i]], col=distinct_hues[5], ci=F);
  lines(fits_ggam[[i]], col=distinct_hues[6], ci=F);
  lines(fits_ggam_orig[[i]], col=distinct_hues[7], ci=F);
  lines(fits_logn[[i]], col=distinct_hues[8], ci=F);
  lines(fits_exp[[i]],col=distinct_hues[9], ci=F);
  lines(fits_genf[[i]],col=distinct_hues[10], ci=F, lty = 2);
  lines(fits_genf_orig[[i]],col=distinct_hues[11], ci=F, lty = 2);  
  lines(fits_rp02[[i]],col=distinct_hues[12], ci=F, lty = 2);
  lines(fits_rp03[[i]],col=distinct_hues[13], ci=F, lty = 2);
  lines(fits_rp04[[i]],col=distinct_hues[14], ci=F, lty = 2);
  lines(fits_rp05[[i]],col=distinct_hues[15], ci=F, lty = 2);
  lines(fits_rp06[[i]],col=distinct_hues[16], ci=F, lty = 2);
  lines(fits_rp07[[i]],col=distinct_hues[17], ci=F, lty = 3);
  lines(fits_rp08[[i]],col=distinct_hues[18], ci=F, lty = 3);
  lines(fits_rp09[[i]],col=distinct_hues[19], ci=F, lty = 3);
  lines(fits_rp10[[i]],col=distinct_hues[20], ci=F, lty = 3)
  
legend("bottomleft", legend = c("Kaplan-Meier","Weibull (AFT)", "Weibull (PH)", "Gompertz", "Log-logistic", "Gamma",
                "Generalized gamma", "Generalized gamma (orig)", "Lognormal", "Exponential", "Generalized F","Generalized F (orig)","RP2","RP3","RP4","RP5", "RP6","RP7", "RP8", "RP9", "RP10"), col = 
         c("red",distinct_hues), 
       title = "Distributions", cex = .95, bty = "n", lty=c(rep(1,10),rep(2,6),rep(3,4)),lwd=3)# lty = 1:2, 
title(main=transition_label[[i]])
}

Figure 19. Vissual Assessment of Survival Curves

endTime <- Sys.time()

paste0("Time in process: ");endTime - tiempo_antes_fits

## [1] "Time in process: "

## Time difference of 3.696909 mins

#23 min aprox.

#For more complicated models, users should specify what covariate values they want summaries for, rather than relying on the default
#</div>
options(warn=0)

if(no_mostrar==1){
jpeg(paste0(dta_path,"_mult_state_ags/exp_surv_int_only.jpg"), height=30, width= 30, res= 600, units = "in")
layout(matrix(1:n_trans, nc = 2, byrow = F))
for (i in 1:n_trans){
plot(km.lc[[i]], col="red", conf.int=F);
  lines(fits_wei[[i]], col=distinct_hues[1], ci=F);
  lines(fits_weiph[[i]], col=distinct_hues[2], ci=F);
  lines(fits_gomp[[i]], col=distinct_hues[3], ci=F);
  lines(fits_llogis[[i]], col=distinct_hues[4], ci=F);#A0A36D
  lines(fits_gam[[i]], col=distinct_hues[5], ci=F);
  lines(fits_ggam[[i]], col=distinct_hues[6], ci=F);
  lines(fits_ggam_orig[[i]], col=distinct_hues[7], ci=F);
  lines(fits_logn[[i]], col=distinct_hues[8], ci=F);
  lines(fits_exp[[i]],col=distinct_hues[9], ci=F);
  lines(fits_genf[[i]],col=distinct_hues[10], ci=F, lty = 2);
  lines(fits_genf_orig[[i]],col=distinct_hues[11], ci=F, lty = 2);  
  lines(fits_rp02[[i]],col=distinct_hues[12], ci=F, lty = 2);
  lines(fits_rp03[[i]],col=distinct_hues[13], ci=F, lty = 2);
  lines(fits_rp04[[i]],col=distinct_hues[14], ci=F, lty = 2);
  lines(fits_rp05[[i]],col=distinct_hues[15], ci=F, lty = 2);
  lines(fits_rp06[[i]],col=distinct_hues[16], ci=F, lty = 2);
  lines(fits_rp07[[i]],col=distinct_hues[17], ci=F, lty = 3);
  lines(fits_rp08[[i]],col=distinct_hues[18], ci=F, lty = 3);
  lines(fits_rp09[[i]],col=distinct_hues[19], ci=F, lty = 3);
  lines(fits_rp10[[i]],col=distinct_hues[20], ci=F, lty = 3)
  
legend("bottomleft", legend = c("Kaplan-Meier","Weibull (AFT)", "Weibull (PH)", "Gompertz", "Log-logistic", "Gamma",
                "Generalized gamma", "Generalized gamma (orig)", "Lognormal", "Exponential", "Generalized F","Generalized F (orig)","RP2","RP3","RP4","RP5", "RP6","RP7", "RP8", "RP9", "RP10"), col = 
         c("red",distinct_hues), 
       title = "Distributions", cex = .95, bty = "n", lty=c(rep(1,10),rep(2,6),rep(3,4)),lwd=3)# lty = 1:2, 
title(main=transition_label[[i]])
}
  dev.off()
}

The following transitions deviated considerably from the models chosen. Consider specify transition parameters manually or check whether to transform the measure of time into an exponential or logarithmic time to event:

• Readmission2 to DWCA3
• Readmission2 to TD3
• Readmission to DWCA2
• Readmission to TD2
• Admission to Readmission
• Admission to DWCA
• Admission to TD

#
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
fit_flexsurvreg0<-data.frame()
fitted_flexsurvreg0<-data.frame()

dists_no_covs_10s<-cbind.data.frame(covs=c(rep("fits_",(20)*n_trans)),
              formal=rep(c("Weibull (AFT)", "Weibull (PH)", "Gompertz", "Log-logistic", "Gamma",
                "Generalized gamma", "Generalized gamma (original)", "Lognormal", "Exponential", "Generalized F", "Generalized F(original)", paste0("RP0",2:9),"RP10"),1*n_trans),
              dist=c("weibull", "weibullph", "llogis", "gamma", "gengamma","gengamma.orig", "gompertz", "lnorm", "exp", "genf","genf.orig",rep("no dist",9)),
              model=rep(c("wei", "weiph", "gomp", "llogis", "gam","ggam", "ggam_orig", "logn", "exp", "genf","genf_orig", paste0("rp0",2:9),"rp10"),1*n_trans),
              trans=rep(1:n_trans, each=20))

for (i in 1:nrow(dists_no_covs_10s)){  #
  
  cat(paste0("#### Flexible Survival Model (w/ covs): ",
               dists_no_covs_10s[i,"formal"], "; transition: ",dists_no_covs_10s[i,"trans"], "\n \n"))  
    
  model<-paste0("fits_",dists_no_covs_10s[i,"model"])
  
  if(!is.null(get(model)[[dists_no_covs_10s[i,"trans"]]])){  
     fitted_flexsurvreg0<- rbind(fitted_flexsurvreg0,cbind.data.frame(dist=rep(dists_no_covs_10s[i,"formal"],), 
                            trans=rep(dists_no_covs_10s[i,"trans"],),
                            data.table::data.table(summary(get(model)[[dists_no_covs_10s[i,"trans"]]], 
                                                           #newdata= newdat2a, 
                                                           #newtime=unique(fitted_km$time), 
                                                           type = "survival", 
                                                           tidy=T)))) 
    # Generate fit indices
    fit_flexsurvreg0<-rbind(fit_flexsurvreg0,
       cbind(dist= dists_no_covs_10s[i,"formal"],
             transition=dists_no_covs_10s[i,"trans"],
             fitstats.flexsurvreg(get(model)[[dists_no_covs_10s[i,"trans"]]])))
    #the BIC may not be appropriate if none of the candidate models are considered to be close to the ‘true’ model.     
    } else {
    cat(paste0("The model that assumed a ",dists_no_covs_10s[i,"formal"]," distribution for the transition number ",dists_no_covs_10s[i,"trans"]," did not converge \n\n"))
    }
}

## #### Flexible Survival Model (w/ covs): Weibull (AFT); transition: 1
##  
## #### Flexible Survival Model (w/ covs): Weibull (PH); transition: 1
##  
## #### Flexible Survival Model (w/ covs): Gompertz; transition: 1
##  
## The model that assumed a Gompertz distribution for the transition number 1 did not converge 
## 
## #### Flexible Survival Model (w/ covs): Log-logistic; transition: 1
##  
## #### Flexible Survival Model (w/ covs): Gamma; transition: 1
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma; transition: 1
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma (original); transition: 1
##  
## #### Flexible Survival Model (w/ covs): Lognormal; transition: 1
##  
## #### Flexible Survival Model (w/ covs): Exponential; transition: 1
##  
## #### Flexible Survival Model (w/ covs): Generalized F; transition: 1
##  
## #### Flexible Survival Model (w/ covs): Generalized F(original); transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP02; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP03; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP04; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP05; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP06; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP07; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP08; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP09; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP10; transition: 1
##  
## #### Flexible Survival Model (w/ covs): Weibull (AFT); transition: 2
##  
## #### Flexible Survival Model (w/ covs): Weibull (PH); transition: 2
##  
## #### Flexible Survival Model (w/ covs): Gompertz; transition: 2
##  
## The model that assumed a Gompertz distribution for the transition number 2 did not converge 
## 
## #### Flexible Survival Model (w/ covs): Log-logistic; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Gamma; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma (original); transition: 2
##  
## #### Flexible Survival Model (w/ covs): Lognormal; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Exponential; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Generalized F; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Generalized F(original); transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP02; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP03; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP04; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP05; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP06; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP07; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP08; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP09; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP10; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Weibull (AFT); transition: 3
##  
## #### Flexible Survival Model (w/ covs): Weibull (PH); transition: 3
##  
## #### Flexible Survival Model (w/ covs): Gompertz; transition: 3
##  
## The model that assumed a Gompertz distribution for the transition number 3 did not converge 
## 
## #### Flexible Survival Model (w/ covs): Log-logistic; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Gamma; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma (original); transition: 3
##  
## #### Flexible Survival Model (w/ covs): Lognormal; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Exponential; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Generalized F; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Generalized F(original); transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP02; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP03; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP04; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP05; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP06; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP07; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP08; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP09; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP10; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Weibull (AFT); transition: 4
##  
## #### Flexible Survival Model (w/ covs): Weibull (PH); transition: 4
##  
## #### Flexible Survival Model (w/ covs): Gompertz; transition: 4
##  
## The model that assumed a Gompertz distribution for the transition number 4 did not converge 
## 
## #### Flexible Survival Model (w/ covs): Log-logistic; transition: 4
##  
## #### Flexible Survival Model (w/ covs): Gamma; transition: 4
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma; transition: 4
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma (original); transition: 4
##  
## #### Flexible Survival Model (w/ covs): Lognormal; transition: 4
##  
## #### Flexible Survival Model (w/ covs): Exponential; transition: 4
##  
## #### Flexible Survival Model (w/ covs): Generalized F; transition: 4
##  
## #### Flexible Survival Model (w/ covs): Generalized F(original); transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP02; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP03; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP04; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP05; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP06; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP07; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP08; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP09; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP10; transition: 4
## 

#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
if(no_mostrar==1){
  fit_flexsurvreg0 %>% 
    dplyr::group_by(transition) %>% 
    summarise(mean(ucl,na.rm=T))
  }
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_

km.lc0<-list()
fitted_km0<-data.frame()

for (i in 1:n_trans){
  km.lc0[[i]] <- survfit(formula= fitform, data = subset(ms_d_match_surv, trans == i))
  fitted_km0<-rbind(fitted_km0,cbind.data.frame(trans=i,fortify(km.lc0[[i]])))
}

#Calculate error
#c("Weibull (AFT)", "Weibull (PH)", "Gompertz", "Log-logistic", "Gamma","Generalized gamma", "Lognormal", "Exponential", "Generalized F")

fitted_flexsurvreg_binned_mix0<-
data.frame(fitted_flexsurvreg0)[,c("dist","trans","time","est","lcl","ucl")] %>% 
  dplyr::left_join(fitted_km0, by=c("trans","time")) %>% 
#there are many observations that was not available due to empirical missing data
#dplyr::filter(!is.na(surv))
  dplyr::group_by(dist,trans) %>% 
  #dplyr::mutate(est= ifelse(is.na(est), mean(est, na.rm=TRUE), est)) %>% 
  #dplyr::mutate(surv= ifelse(is.na(surv), mean(surv, na.rm=TRUE), surv)) %>% 
  tidyr::fill(surv,.direction="updown") %>% 
  tidyr::fill(est,.direction="updown") %>% 
  ungroup()

db_for_apply_rmse0<-
  cbind.data.frame(
      dist=rep(c("Weibull (AFT)", "Weibull (PH)", "Gompertz", "Log-logistic", "Gamma",
                "Generalized gamma", "Generalized gamma (original)", "Lognormal", "Exponential", "Generalized F", "Generalized F (orignial)", paste0("RP0",2:9),"RP10"),n_trans),
      trans=rep(c(1:n_trans),each=20*2))
   
rmse_comp_fits0<- data.frame()
for(i in 1:nrow(db_for_apply_rmse0)){
  rmse<- Metrics::rmse(subset(fitted_flexsurvreg_binned_mix0, dist==db_for_apply_rmse0[i,"dist"] & 
                       trans==db_for_apply_rmse0[i,"trans"])$est,
              subset(fitted_flexsurvreg_binned_mix0, dist==db_for_apply_rmse0[i,"dist"] & 
                       trans==db_for_apply_rmse0[i,"trans"])$surv)

  rmse_comp_fits0<- rbind(rmse_comp_fits0,cbind(dist=db_for_apply_rmse0[i,"dist"],
                                                  residential=db_for_apply_rmse0[i,"residential"],
                                                  trans=db_for_apply_rmse0[i,"trans"],
                                                  rmse=rmse))
}

rmse_comp_fits0_filter<-
  rmse_comp_fits0 %>% 
    dplyr::filter(rmse!="NaN") %>%  
      dplyr::arrange(trans,dist, rmse)%>%
      dplyr::mutate(rmse=round(as.numeric(rmse),4))

setting_type_label<- c(`0`="Outpatient",`1`="Residential") 
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_

We checked the AICs, cumulative relative error (differences between the Kaplan-Meier curve, and the observed behavior of each distribution).

#load("C:/Users/CISS Fondecyt/OneDrive/Escritorio/mult_state_ags.RData")

c26 <- c(
  "dodgerblue2", # red
  "green4",
  "#6A3D9A", # purple
  "#FF7F00", # orange
  "gray16", "gold1",
  "skyblue2", "#FB9A99", # lt pink
  "palegreen2",
  "#CAB2D6", # lt purple
  "#FDBF6F", # lt orange
  "gray70", "khaki2",
  "maroon", "orchid1", "deeppink1", "blue1", "steelblue4",
  "darkturquoise", "green1", "yellow4", "yellow3",
  "darkorange4", "brown", "gray40")

fitted_flexsurvreg_binned_mix0_plot<-
fitted_flexsurvreg_binned_mix0 %>% 
  dplyr::mutate(abs((est-surv)/surv)) %>% group_by(dist,trans) %>% dplyr::mutate(cumsum_error=cumsum(`abs((est - surv)/surv)`)) %>% 
  dplyr::ungroup() %>% 
  dplyr::mutate(text=paste0("Distribution: ",dist,"<br>survival: ",sprintf("%1.2f",est),"<br>Time (in years): ",round(time/365.25,2), "<br>CIs: ",ifelse(!is.na(lcl),"CIs","no CIs"),"<br>Cumsum(abs((est-surv)/surv))= ",round(cumsum_error,2))) %>% 
  dplyr::mutate(text2=paste0("Distribution: KM","<br>survival: ",sprintf("%1.2f",surv),"<br>Time (in years): ",round(time/365.25,2))) %>%
ggplot()+
    geom_line(aes(time, est, color=dist),size=1)+
    geom_line(aes(time, surv), color="red",size=1)+
    scale_color_manual(name="Distributions", values = c26[1:20])+
                         #c("black",brewer.pal(n = 9, name = 'Paired')))+
                         #c("#112A60","#085754","#D3A347","#4F3C91","red","#112A60","#085754","#8F630D","#251363")) +
    facet_wrap(.~trans,labeller = labeller(trans = transition_label), ncol=3, scales = "free_y")+
    sjPlot::theme_sjplot2()+
    theme(legend.position="bottom",
          strip.background = element_rect(fill = "white", colour = "white"))+
  scale_x_continuous(breaks = seq(0, 365.25*12, 2*365.25), labels=seq(0,12,2))+
  #ylim(0,1.25)+
  #theme(axis.text.x = element_blank(), 
  #      panel.grid.major = element_blank(), 
  #      panel.grid.minor = element_blank()) +
  labs(y="",x="")+
  theme(legend.position='none')

m <- list(
  l = 80*1.05,
  r = 80,
  b = 80,
  t = 80*1.05,
  pad = 4
)
p1<-
fitted_flexsurvreg_binned_mix0_plot$data %>% 
  dplyr::mutate(years=round(time/365.25,0)) %>% 
  dplyr::mutate(trans2=trans) %>% 
  dplyr::left_join(cbind.data.frame(transition_label,trans2=1:n_trans), by="trans2") %>% 
group_by(trans) %>%
group_map(~ plot_ly(data=., x = ~time, y = ~est, color = ~dist, type = "scatter", mode="lines", text=~text) %>% 
            add_trace(y = ~surv, type = 'scatter', mode = 'lines', line = list(color = 'red', width = 1), text=~text2)%>% 
  layout(annotations = list(list(x = 0.5 , y = 1.03, text = ~unique(transition_label), showarrow = F, xref='paper', yref='paper'))) %>% 
  layout(
    xaxis = list(
      ticktext = list(1, 3, 5, 7, 9, 11), 
      textangle = 0,
      tickvals = list(365.25*1, 365.25*3, 365.25*5, 365.25*7, 365.25*9, 365.25*11),
      tickmode = "array"
  )), keep=TRUE) %>%    
subplot(nrows = 3, shareX = T, shareY=T, titleX = F, titleY = F, margin = .05)%>% layout(showlegend = F) %>% #, margin = 0.05) %>% 
  layout(annotations = list(
                list(x = -0.03, y = 0.5, text = "Survival",
                     font = list(color = "black",size = 15),
                     textangle = 270,
                     showarrow = F, xref='paper', yref='paper', size=48,margin =m)))

p1

Figure 20. Vissual Assessment of Survival Curves (w/o covars)

  #htmlwidgets::saveWidget(as_widget(partial_bundle(p1, local=T)), "test1.html")

#https://github.com/ropensci/plotly/issues/1586

                 #%>% 
  #layout(annotations = list(
  #              list(x = 0.5 , y = -0.03, text = "Survival",
  #                   font = list(color = "black",size = 15),
  #                   textangle = 0,
  #                   showarrow = F, xref='paper', yref='paper', size=48,
  #                   automargin = T)
   #             ))

n_trans2 <- max(trans_matrix, na.rm = T)

fit_flexurvreg0_kable<-
fit_flexsurvreg0 %>% 
  dplyr::arrange(transition, AIC) %>% 
  dplyr::left_join(cbind.data.frame(transition_label,trans_nmb=1:n_trans2),by=c("transition"="trans_nmb")) %>% 
  dplyr::mutate(trans_w_nmb=paste0(sprintf("%02d",transition),")",transition_label)) %>% 
  dplyr::select(trans_w_nmb, dist, Df, n2ll, AIC, AICc, BIC) %>%
  dplyr::mutate_at(vars(n2ll, AIC, AICc, BIC),~ifelse(abs(.)>1e6,NA,format(round(.,2), scientific = FALSE)))
  
 
fit_flexurvreg0_kable %>%   
      knitr::kable(.,format = "html", format.args = list(decimal.mark = ".", big.mark = ","),
               caption = paste0("Table 11. Fit indices of the adjusted survival analyses (intercept-only)"),
               col.names = c("Transition","Distribution", "DF","Negative 2 Log Likelihood","AIC","AICc","BIC"),
               align =c("l",rep('c', 101))) %>%
  kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12) %>%
  #kableExtra::add_footnote("Note. NA= Null values", notation="none") %>% 
  kableExtra::scroll_box(width = "100%", height = "375px") 

Table 11. Fit indices of the adjusted survival analyses (intercept-only)

Transition	Distribution	DF	Negative 2 Log Likelihood	AIC	AICc	BIC
01)Admission to Readmission	Generalized F(original)	4	NA	NA	NA	NA
01)Admission to Readmission	RP06	8	121036.84	121052.84	121052.85	121117.16
01)Admission to Readmission	RP04	6	121041.85	121053.85	121053.86	121102.09
01)Admission to Readmission	RP08	10	121033.99	121053.99	121054.00	121134.39
01)Admission to Readmission	Generalized F	4	121046.51	121054.51	121054.51	121086.67
01)Admission to Readmission	RP10	12	121031.52	121055.52	121055.53	121151.99
01)Admission to Readmission	RP09	11	121033.96	121055.96	121055.97	121144.39
01)Admission to Readmission	RP07	9	121038.72	121056.72	121056.73	121129.08
01)Admission to Readmission	RP05	7	121043.64	121057.64	121057.65	121113.92
01)Admission to Readmission	RP03	5	121058.64	121068.64	121068.65	121108.84
01)Admission to Readmission	RP02	4	121068.72	121076.72	121076.72	121108.88
01)Admission to Readmission	Generalized gamma	3	121291.72	121297.72	121297.72	121321.84
01)Admission to Readmission	Lognormal	2	121619.49	121623.49	121623.49	121639.57
01)Admission to Readmission	Generalized gamma (original)	3	121759.04	121765.04	121765.05	121789.16
01)Admission to Readmission	Log-logistic	2	122122.85	122126.85	122126.85	122142.93
01)Admission to Readmission	Weibull (AFT)	2	122464.79	122468.79	122468.79	122484.87
01)Admission to Readmission	Weibull (PH)	2	122464.79	122468.79	122468.79	122484.87
01)Admission to Readmission	Gamma	2	122578.49	122582.49	122582.49	122598.57
01)Admission to Readmission	Exponential	1	122964.89	122966.89	122966.89	122974.93
02)Readmission to Readmission2	Generalized F	4	36551.44	36559.44	36559.45	36586.50
02)Readmission to Readmission2	Generalized F(original)	4	36551.44	36559.44	36559.45	36586.50
02)Readmission to Readmission2	RP03	5	36551.19	36561.19	36561.20	36595.01
02)Readmission to Readmission2	RP04	6	36549.46	36561.46	36561.48	36602.05
02)Readmission to Readmission2	RP05	7	36547.95	36561.95	36561.96	36609.29
02)Readmission to Readmission2	RP08	10	36543.34	36563.34	36563.37	36630.97
02)Readmission to Readmission2	RP09	11	36541.54	36563.54	36563.58	36637.94
02)Readmission to Readmission2	RP06	8	36547.72	36563.72	36563.74	36617.83
02)Readmission to Readmission2	RP07	9	36545.91	36563.91	36563.94	36624.78
02)Readmission to Readmission2	RP10	12	36541.57	36565.57	36565.62	36646.73
02)Readmission to Readmission2	RP02	4	36561.33	36569.33	36569.34	36596.39
02)Readmission to Readmission2	Generalized gamma	3	36651.53	36657.53	36657.54	36677.82
02)Readmission to Readmission2	Lognormal	2	36676.37	36680.37	36680.37	36693.90
02)Readmission to Readmission2	Generalized gamma (original)	3	36704.13	36710.13	36710.13	36730.42
02)Readmission to Readmission2	Log-logistic	2	36786.43	36790.43	36790.43	36803.96
02)Readmission to Readmission2	Weibull (AFT)	2	36905.75	36909.75	36909.76	36923.28
02)Readmission to Readmission2	Weibull (PH)	2	36905.75	36909.75	36909.76	36923.28
02)Readmission to Readmission2	Gamma	2	36931.00	36935.00	36935.00	36948.53
02)Readmission to Readmission2	Exponential	1	36978.75	36980.75	36980.75	36987.51
03)Readmission2 to Readmission3	Generalized F(original)	4	NA	NA	NA	NA
03)Readmission2 to Readmission3	Generalized F	4	11510.25	11518.25	11518.27	11540.68
03)Readmission2 to Readmission3	RP02	4	11511.22	11519.22	11519.24	11541.64
03)Readmission2 to Readmission3	RP03	5	11511.36	11521.36	11521.39	11549.39
03)Readmission2 to Readmission3	RP04	6	11510.57	11522.57	11522.61	11556.21
03)Readmission2 to Readmission3	RP05	7	11510.45	11524.45	11524.51	11563.70
03)Readmission2 to Readmission3	RP06	8	11509.57	11525.57	11525.64	11570.42
03)Readmission2 to Readmission3	RP07	9	11509.11	11527.11	11527.20	11577.56
03)Readmission2 to Readmission3	RP08	10	11507.51	11527.51	11527.62	11583.57
03)Readmission2 to Readmission3	RP09	11	11506.53	11528.53	11528.66	11590.20
03)Readmission2 to Readmission3	RP10	12	11505.30	11529.30	11529.46	11596.58
03)Readmission2 to Readmission3	Generalized gamma	3	11532.03	11538.03	11538.04	11554.85
03)Readmission2 to Readmission3	Lognormal	2	11541.52	11545.52	11545.53	11556.73
03)Readmission2 to Readmission3	Generalized gamma (original)	3	11552.33	11558.33	11558.34	11575.15
03)Readmission2 to Readmission3	Log-logistic	2	11578.52	11582.52	11582.53	11593.73
03)Readmission2 to Readmission3	Weibull (AFT)	2	11621.91	11625.91	11625.92	11637.12
03)Readmission2 to Readmission3	Weibull (PH)	2	11621.91	11625.91	11625.92	11637.12
03)Readmission2 to Readmission3	Gamma	2	11626.91	11630.91	11630.91	11642.12
03)Readmission2 to Readmission3	Exponential	1	11630.65	11632.65	11632.66	11638.26
04)Readmission3 to Readmission4	Generalized F(original)	4	3607.60	3615.60	3615.67	3633.49
04)Readmission3 to Readmission4	Generalized F	4	3607.67	3615.67	3615.73	3633.56
04)Readmission3 to Readmission4	RP02	4	3608.19	3616.19	3616.25	3634.08
04)Readmission3 to Readmission4	RP03	5	3608.33	3618.33	3618.42	3640.69
04)Readmission3 to Readmission4	RP04	6	3607.92	3619.92	3620.05	3646.75
04)Readmission3 to Readmission4	RP05	7	3606.91	3620.91	3621.08	3652.21
04)Readmission3 to Readmission4	RP06	8	3606.40	3622.40	3622.63	3658.18
04)Readmission3 to Readmission4	Lognormal	2	3619.93	3623.93	3623.95	3632.88
04)Readmission3 to Readmission4	RP07	9	3606.32	3624.32	3624.60	3664.57
04)Readmission3 to Readmission4	Generalized gamma	3	3619.45	3625.45	3625.48	3638.86
04)Readmission3 to Readmission4	RP08	10	3605.47	3625.47	3625.82	3670.20
04)Readmission3 to Readmission4	RP09	11	3605.77	3627.77	3628.18	3676.96
04)Readmission3 to Readmission4	Generalized gamma (original)	3	3622.23	3628.23	3628.27	3641.65
04)Readmission3 to Readmission4	RP10	12	3605.95	3629.95	3630.44	3683.61
04)Readmission3 to Readmission4	Log-logistic	2	3627.00	3631.00	3631.02	3639.95
04)Readmission3 to Readmission4	Exponential	1	3641.53	3643.53	3643.54	3648.00
04)Readmission3 to Readmission4	Weibull (AFT)	2	3640.15	3644.15	3644.17	3653.10
04)Readmission3 to Readmission4	Weibull (PH)	2	3640.15	3644.15	3644.17	3653.10
04)Readmission3 to Readmission4	Gamma	2	3641.13	3645.13	3645.14	3654.07

rmse_plot0<-
rmse_comp_fits0_filter %>% 
  dplyr::mutate(trans=as.numeric(trans)) %>% 
  dplyr::mutate(trans2=as.numeric(trans)) %>% 
    dplyr::left_join(cbind.data.frame(transition_label,trans_nmb=1:n_trans2),by=c("trans"="trans_nmb")) %>% 
  dplyr::mutate(trans_w_nmb=paste0(sprintf("%02d",trans),")",transition_label)) %>% 
        
   ggplot()+
  geom_bar(aes(x=dist, y=rmse), position="dodge", stat="identity", alpha=0.4)+
  sjPlot::theme_sjplot2()+
  #geom_errorbar(aes(x=t, ymin=L, ymax=U, color=program), position="dodge", stat="identity", width=0.4,  alpha=0.9, size=1.3)+
  facet_wrap(.~trans_w_nmb, ncol=3, scales="free_y", dir="h") + 
  xlab("") + 
  ylab("")+
  #ylab("State occupancy probabilities") + 
  theme_minimal()+
  theme(legend.position="bottom")+
    theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        strip.background = element_blank(),
        panel.border = element_blank(),
          legend.box.background = element_blank())+
  theme(axis.text.x = element_text(angle = 90,hjust=0.95,vjust=0.2))#, vjust=0, hjust=0.3

rmse_plot0

Figure 21. Suvival RMSEs, Ten-states Model (w/o covars)

#install.packages("survHE")
library(survHE)
#First, defines the vector of models to be used
mods <- c("weibull", "weibullph", "llogis", "gamma", "gengamma", "gompertz", "lnorm" ,"exp")#, "genf")

# And then runs the models using MLE via flexsurv
#m2 <- fit.models(formula = formula, data = data, distr = c("exp","wei","lno","llo"), method="inla")
#m3 <- fit.models(formula = formula, data = subset(ms_d_match_surv, trans == 1), distr = mods, method="hmc")
#Error in model.matrix(formula, data)[(mf %>% filter(event == 1))$ID, ] : 
#  subíndice fuera de  los límites

list_aics<-list()
for (i in 1:n_trans){
  m1 <- fit.models(formula = fitform, data = subset(ms_d_match_surv, trans == i), distr = mods)
  model.fit.plot(MLE=m1, stacked = T, #HMC=m3, 
               models = c("Weibull (AFT)", "Weibull (PH)", "Gompertz", "Log-logistic", "Gamma",
                "Generalized gamma", "Lognormal", "Exponential", "Generalized F"),
               #lab.profile = c("Outpatient","Residential")
               name_legend = "Inferential method")+ 
  ggtitle("Model comparison based on AIC")+
  labs(subtitle=plots[i,"title"])
}

tiempo_antes_fits2<-Sys.time()

newtime0 = seq(from=round(min(ms_d_match_surv$time),0), to= round(max(ms_d_match_surv$time),0))

#_#_#_#_#_#_#_#_#_#_#_#_#_
#covariates

#Specify the formula
fitform2 <- Surv(time, status) ~  factor(tipo_de_plan_res_1)


kernel_haz_est2a<-list()
kernel_haz_est2b<-list()
kernel_haz2a<-list()
kernel_haz2b<-list()
for (i in 1:n_trans){
library("muhaz")
kernel_haz_est2a[[i]] <- muhaz(ms_d_match_surv[which(ms_d_match_surv$trans==i &
                        ms_d_match_surv$tipo_de_plan_res_1==1),"time"],
                        ms_d_match_surv[which(ms_d_match_surv$trans==i &
                        ms_d_match_surv$tipo_de_plan_res_1==1),"status"])
kernel_haz2a[[i]] <- data.table(time = kernel_haz_est2a[[i]]$est.grid,
                         est = kernel_haz_est2a[[i]]$haz.est,
                         dist = "Kernel density")
kernel_haz_est2b[[i]] <- muhaz(ms_d_match_surv[which(ms_d_match_surv$trans==i &
                        ms_d_match_surv$tipo_de_plan_res_1==0),"time"],
                        ms_d_match_surv[which(ms_d_match_surv$trans==i &
                        ms_d_match_surv$tipo_de_plan_res_1==0),"status"])
kernel_haz2b[[i]] <- data.table(time = kernel_haz_est2b[[i]]$est.grid,
                         est = kernel_haz_est2b[[i]]$haz.est,
                         dist = "Kernel density")
}

haz_int_only2<-
  rbind(cbind(trans=rep(1:n_trans,each=nrow(kernel_haz2a[[i]])),
              tipo_de_plan_res_1=rep(1,nrow(kernel_haz2a[[i]])),
      rbindlist(kernel_haz2a)),
      cbind(trans=rep(1:n_trans,each=nrow(kernel_haz2b[[i]])),
            tipo_de_plan_res_1=rep(0,nrow(kernel_haz2b[[i]])),
      rbindlist(kernel_haz2b)))

#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_

n_trans2 <- max(trans_matrix, na.rm = T)

dists_w_covs_10s<-cbind.data.frame(covs=c(rep("fits_c_",(20)*n_trans2)),
              formal=rep(c("Weibull (AFT)", "Weibull (PH)", "Gompertz", "Log-logistic", "Gamma",
                "Generalized gamma", "Generalized gamma (original)", "Lognormal", "Exponential", "Generalized F", "Generalized F(original)", paste0("RP0",2:9),"RP10"),1*n_trans2),
              dist=c("weibull", "weibullph", "llogis", "gamma", "gengamma","gengamma.orig", "gompertz", "lnorm", "exp", "genf","genf.orig",rep("no dist",9)),
              model=rep(c("wei", "weiph", "gomp", "llogis", "gam","ggam", "ggam_orig", "logn", "exp", "genf","genf_orig", paste0("rp0",2:9),"rp10"),1*n_trans2),
              trans=rep(1:n_trans2, each=20))

#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
no_attempts <- 30

km.lc2a<-list()
km.lc2b<-list()
fits_cox2<-list()
fits_c_cox2<-list()

for (i in 1:n_trans2){
  r <- NULL
  attempt <- 0
  while( is.null(r) && attempt <= no_attempts ) {
    attempt <- attempt + 1
    try(
      r <- flexsurvreg(formula=fitform2,
                                 data = subset(ms_d_match_surv, trans == i),
                                 dist = "weibull")
    )
  } 
  fits_c_wei[[i]] <- r
}

for (i in 1:n_trans2){
  r <- NULL
  attempt <- 0
  while( is.null(r) && attempt <= no_attempts ) {
    attempt <- attempt + 1
    try(
      r <- flexsurvreg(formula=fitform2,
                                 data = subset(ms_d_match_surv, trans == i),
                                 dist = "weibullph")
    )
  } 
  fits_c_weiph[[i]] <- r
}

for (i in 1:n_trans2){
  r <- NULL
  attempt <- 0
  while( is.null(r) && attempt <= no_attempts ) {
    attempt <- attempt + 1
    try(
      r <- flexsurvreg(formula=fitform2,
                                 data = subset(ms_d_match_surv, trans == i),
                                 dist = "llogis")
    )
  } 
  fits_c_llogis[[i]] <- r
}

for (i in 1:n_trans2){
  r <- NULL
  attempt <- 0
  while( is.null(r) && attempt <= no_attempts ) {
    attempt <- attempt + 1
    try(
      r <- flexsurvreg(formula=fitform2,
                                 data = subset(ms_d_match_surv, trans == i),
                                 dist = "gamma")
    )
  } 
  fits_c_gam[[i]] <- r
}

for (i in 1:n_trans2){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvreg(formula=fitform2,
                                   data = subset(ms_d_match_surv, trans == i),
                                   dist = "gengamma")
      )
    }
    fits_c_ggam[[i]] <- r
}

for (i in 1:n_trans2){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvreg(formula=fitform2,
                                   data = subset(ms_d_match_surv, trans == i), 
                                   dist = "gompertz")
      )
    }
    fits_c_gomp[[i]] <- r
}  

## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito

for (i in 1:n_trans2){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvreg(formula=fitform2,
                                   data = subset(ms_d_match_surv, trans == i),
                                   dist = "lnorm")
      )
    }
    fits_c_logn[[i]] <- r
}  

for (i in 1:n_trans2){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvreg(formula=fitform2,
                                   data = subset(ms_d_match_surv, trans == i),
                                   dist = "exp")
      )
    }
    fits_c_exp[[i]] <- r
}  


for (i in 1:n_trans2){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvreg(formula=fitform2,
                                   data = subset(ms_d_match_surv, trans == i),
                                   dist = "genf")
      )
    }
    fits_c_genf[[i]] <- r
}  


for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
        attempt <- attempt + 1
        try(
            r <- flexsurvreg(formula=fitform2,
                             data = subset(ms_d_match_surv, trans == i),
                             dist = "gengamma.orig")
        )
    }
    fits_c_ggam_orig[[i]] <- r
}  

for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
        attempt <- attempt + 1
        try(
            r <- flexsurvreg(formula=fitform2,
                             data = subset(ms_d_match_surv, trans == i),
                             dist = "genf.orig")
        )
    }
    fits_c_genf_orig[[i]] <- r
}

## Warning in flexsurvreg(formula = fitform2, data = subset(ms_d_match_surv, :
## Optimisation has probably not converged to the maximum likelihood - Hessian is
## not positive definite.

## Warning in flexsurvreg(formula = fitform2, data = subset(ms_d_match_surv, :
## Optimisation has probably not converged to the maximum likelihood - Hessian is
## not positive definite.

for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform2,k=2,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_c_rp02[[i]] <- r
}  
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform2,k=3,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_c_rp03[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform2,k=4,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_c_rp04[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform2,k=5,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_c_rp05[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform2,k=6,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_c_rp06[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform2,k=7,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_c_rp07[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform2,k=8,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_c_rp08[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform2,k=9,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_c_rp09[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform2,k=10,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_c_rp10[[i]] <- r
}

#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_

fitted_km<-data.frame()
for (i in 1:n_trans2){
km.lc2a[[i]] <- survfit(formula= fitform, data = subset(ms_d_match_surv, trans == i & tipo_de_plan_res_1==1))
km.lc2b[[i]] <- survfit(formula= fitform, data = subset(ms_d_match_surv, trans == i & tipo_de_plan_res_1==0))

fitted_km<-rbind(fitted_km,cbind.data.frame(trans=i,residential=rep(1,), fortify(km.lc2a[[i]])))
fitted_km<-rbind(fitted_km,cbind.data.frame(trans=i,residential=rep(0,), fortify(km.lc2b[[i]])))
}

#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_

newdat2a <- data.table::data.table(tipo_de_plan_res_1= factor(c(rep(1,1))))
newdat2b <- data.table::data.table(tipo_de_plan_res_1= factor(c(rep(0,1))))

fitted_flexsurvreg<-data.frame()
fit_flexsurvreg<-data.frame()

for (i in 1:nrow(dists_w_covs_10s)){  #
  
cat(paste0("#### Flexible Survival Model (w/ covs): ",
             dists_w_covs_10s[i,"formal"], "; transition: ",dists_w_covs_10s[i,"trans"], "\n \n"))  
  
model<-paste0("fits_c_",dists_w_covs_10s[i,"model"])

if(!is.null(get(model)[[dists_w_covs_10s[i,"trans"]]])){  
  #Generate databases
 fitted_flexsurvreg<- rbind(fitted_flexsurvreg,cbind.data.frame(dist=rep(dists_w_covs_10s[i,"formal"],), 
                            trans=rep(dists_w_covs_10s[i,"trans"],),
                            residential=rep(1,),
                            data.table::data.table(summary(get(model)[[dists_w_covs_10s[i,"trans"]]], newdata= newdat2a, newtime=unique(fitted_km$time), type = "survival", tidy=T)))) 
  fitted_flexsurvreg<- rbind(fitted_flexsurvreg,cbind.data.frame(dist=rep(dists_w_covs_10s[i,"formal"],), 
                            trans=rep(dists_w_covs_10s[i,"trans"],),
                            residential=rep(0,),
                            data.table::data.table(summary(get(model)[[dists_w_covs_10s[i,"trans"]]], newdata= newdat2b,  newtime=unique(fitted_km$time), type = "survival", tidy=T)))) 
 #t=newtime0, 
 
  # Generate fit indices
  fit_flexsurvreg<-rbind(fit_flexsurvreg,
     cbind(dist= dists_w_covs_10s[i,"formal"],
           transition=dists_w_covs_10s[i,"trans"],
           fitstats.flexsurvreg(get(model)[[dists_w_covs_10s[i,"trans"]]])))
  #the BIC may not be appropriate if none of the candidate models are considered to be close to the ‘true’ model.     
  } else {
  cat(paste0("The model that assumed a ",dists_w_covs_10s[i,"formal"]," distribution for the transition number ",dists_w_covs_10s[i,"trans"]," did not converge \n\n"))
  }
}

## #### Flexible Survival Model (w/ covs): Weibull (AFT); transition: 1
##  
## #### Flexible Survival Model (w/ covs): Weibull (PH); transition: 1
##  
## #### Flexible Survival Model (w/ covs): Gompertz; transition: 1
##  
## The model that assumed a Gompertz distribution for the transition number 1 did not converge 
## 
## #### Flexible Survival Model (w/ covs): Log-logistic; transition: 1
##  
## #### Flexible Survival Model (w/ covs): Gamma; transition: 1
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma; transition: 1
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma (original); transition: 1
##  
## #### Flexible Survival Model (w/ covs): Lognormal; transition: 1
##  
## #### Flexible Survival Model (w/ covs): Exponential; transition: 1
##  
## #### Flexible Survival Model (w/ covs): Generalized F; transition: 1
##  
## #### Flexible Survival Model (w/ covs): Generalized F(original); transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP02; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP03; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP04; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP05; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP06; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP07; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP08; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP09; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP10; transition: 1
##  
## #### Flexible Survival Model (w/ covs): Weibull (AFT); transition: 2
##  
## #### Flexible Survival Model (w/ covs): Weibull (PH); transition: 2
##  
## #### Flexible Survival Model (w/ covs): Gompertz; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Log-logistic; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Gamma; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma (original); transition: 2
##  
## #### Flexible Survival Model (w/ covs): Lognormal; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Exponential; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Generalized F; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Generalized F(original); transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP02; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP03; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP04; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP05; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP06; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP07; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP08; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP09; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP10; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Weibull (AFT); transition: 3
##  
## #### Flexible Survival Model (w/ covs): Weibull (PH); transition: 3
##  
## #### Flexible Survival Model (w/ covs): Gompertz; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Log-logistic; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Gamma; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma (original); transition: 3
##  
## #### Flexible Survival Model (w/ covs): Lognormal; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Exponential; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Generalized F; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Generalized F(original); transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP02; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP03; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP04; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP05; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP06; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP07; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP08; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP09; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP10; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Weibull (AFT); transition: 4
##  
## #### Flexible Survival Model (w/ covs): Weibull (PH); transition: 4
##  
## #### Flexible Survival Model (w/ covs): Gompertz; transition: 4
##  
## #### Flexible Survival Model (w/ covs): Log-logistic; transition: 4
##  
## #### Flexible Survival Model (w/ covs): Gamma; transition: 4
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma; transition: 4
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma (original); transition: 4
##  
## #### Flexible Survival Model (w/ covs): Lognormal; transition: 4
##  
## #### Flexible Survival Model (w/ covs): Exponential; transition: 4
##  
## #### Flexible Survival Model (w/ covs): Generalized F; transition: 4
##  
## #### Flexible Survival Model (w/ covs): Generalized F(original); transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP02; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP03; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP04; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP05; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP06; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP07; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP08; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP09; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP10; transition: 4
## 

#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
if(no_mostrar==1){
  fit_flexsurvreg %>% 
    dplyr::group_by(trans) %>% 
    summarise(mean(ucl,na.rm=T))
  }

#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_

#load("C:/Users/CISS Fondecyt/OneDrive/Escritorio/mult_state_ags.RData")

c26 <- c(
  "dodgerblue2", # red
  "green4",
  "#6A3D9A", # purple
  "#FF7F00", # orange
  "gray16", "gold1",
  "skyblue2", "#FB9A99", # lt pink
  "palegreen2",
  "#CAB2D6", # lt purple
  "#FDBF6F", # lt orange
  "gray70", "khaki2",
  "maroon", "orchid1", "deeppink1", "blue1", "steelblue4",
  "darkturquoise", "green1", "yellow4", "yellow3",
  "darkorange4", "brown", "gray40")

fitted_flexsurvreg_plot<-
fitted_flexsurvreg %>%  
  dplyr::left_join(fitted_km, by=c("time", "trans", "residential")) %>% 
  group_by(trans,dist,residential) %>% 
  tidyr::fill(surv,.direction="down") %>% 
  ungroup() %>% 
  dplyr::mutate(abs((est-surv)/surv)) %>% 
  group_by(trans,dist,residential) %>% 
  dplyr::mutate(cumsum_error=cumsum(replace_na(`abs((est - surv)/surv)`, 0))) %>% 
  dplyr::ungroup() %>% 
  dplyr::mutate(text=paste0("Distribution: ",dist,"<br>survival: ",sprintf("%1.2f",est),"<br>Time (in years): ",round(time/365.25,2), "<br>CIs: ",ifelse(!is.na(lcl),"CIs","no CIs"),"<br>Cumsum(abs((est-surv)/surv))= ",round(cumsum_error,2),"<br>Residential: ",ifelse(residential==1,"Yes","No"))) %>% 
  dplyr::mutate(text2=paste0("Distribution: KM","<br>survival: ",sprintf("%1.2f",surv),"<br>Time (in years): ",round(time/365.25,2),"<br>Residential: ",ifelse(residential==1,"Yes","No"))) %>%
ggplot()+
    geom_line(aes(time, est, color=dist, linetype=factor(residential)),size=1)+
    geom_line(aes(time, surv, linetype=factor(residential)), color="red",size=1)+
    scale_color_manual(name="Distributions", values = c26[1:20])+
    facet_wrap(.~trans,labeller = labeller(trans = transition_label), ncol=3, scales = "free_y")+
    sjPlot::theme_sjplot2()+
    theme(legend.position="bottom",
          strip.background = element_rect(fill = "white", colour = "white"))+
 # scale_x_continuous(breaks = seq(0, 365.25*12, 2*365.25), labels=seq(0,12,2))+
  labs(y="",x="")+
  theme(legend.position='none')

m <- list(
  l = 80*1.05,
  r = 80,
  b = 80,
  t = 80*1.05,
  pad = 4
)

p2<-
fitted_flexsurvreg_plot$data %>% 
  dplyr::mutate(years=round(time/365.25,0)) %>% 
  dplyr::mutate(trans2=trans) %>% 
  dplyr::left_join(cbind.data.frame(transition_label,trans2=1:n_trans2), by="trans2") %>% 
group_by(trans) %>%
group_map(~ plot_ly(data=., x = ~time, y = ~est, color = ~dist, type = "scatter", mode="lines", linetype = ~residential, text=~text) %>% 
            add_trace(y = ~surv, type = 'scatter', mode = 'lines', line = list(color = 'red', width = 1, linetype = ~residential), text=~text2)%>% 
  layout(annotations = list(list(x = 0.5 , y = 1.03, text = ~unique(transition_label), showarrow = F, xref='paper', yref='paper'))) %>% 
  layout(
    xaxis = list(
      ticktext = list(1, 3, 5, 7, 9, 11), 
      textangle = 0,
      tickvals = list(365.25*1, 365.25*3, 365.25*5, 365.25*7, 365.25*9, 365.25*11),
      tickmode = "array"
  )), keep=TRUE) %>%    
subplot(nrows = 3, shareX = T, shareY=T, titleX = F, titleY = F, margin = .05)%>% layout(showlegend = F) %>% #, margin = 0.05) %>% 
  layout(annotations = list(
                list(x = -0.1, y = 0.5, text = "Survival",
                     font = list(color = "black",size = 15),
                     textangle = 270,
                     showarrow = F, xref='paper', yref='paper', size=48,margin =m)))

## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors

## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors

## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors

## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors

## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors

## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors

## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors

## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors

#p2

  htmlwidgets::saveWidget(as_widget(p2), "test2_jun.html")
#https://docs.ropensci.org/plotly/reference/partial_bundle.html

File is available in this link.

fit_flexurvreg_kable<-
fit_flexsurvreg %>% 
  dplyr::arrange(transition, AIC) %>% 
  dplyr::left_join(cbind.data.frame(transition_label,trans_nmb=1:n_trans2),by=c("transition"="trans_nmb")) %>% 
  dplyr::mutate(trans_w_nmb=paste0(sprintf("%02d",transition),")",transition_label)) %>% 
  dplyr::select(trans_w_nmb,dist,Df,n2ll,AIC,AICc,BIC)%>%
  dplyr::mutate_at(vars(n2ll, AIC, AICc, BIC),~ifelse(abs(.)>1e6,NA,format(round(.,2), scientific = FALSE)))  
  
fit_flexurvreg_kable %>%   
      knitr::kable(.,format = "html", format.args = list(decimal.mark = ".", big.mark = ","),
               caption = paste0("Table 12. Fit indices of the adjusted survival analyses"),
               col.names = c("Transition","Distribution", "DF","Negative 2 Log Likelihood","AIC","AICc","BIC"),
               align =c("l",rep('c', 101))) %>%
  kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12) %>%
  #kableExtra::add_footnote("Note. NA= Null values", notation="none") %>% 
  kableExtra::scroll_box(width = "100%", height = "375px") 

Table 12. Fit indices of the adjusted survival analyses

Transition	Distribution	DF	Negative 2 Log Likelihood	AIC	AICc	BIC
01)Admission to Readmission	Generalized F(original)	5	NA	NA	NA	NA
01)Admission to Readmission	Generalized F	5	120769.03	120779.03	120779.03	120819.23
01)Admission to Readmission	RP06	9	120779.16	120797.16	120797.17	120869.52
01)Admission to Readmission	RP04	7	120784.24	120798.24	120798.25	120854.52
01)Admission to Readmission	RP08	11	120776.37	120798.37	120798.38	120886.81
01)Admission to Readmission	RP10	13	120773.93	120799.93	120799.95	120904.45
01)Admission to Readmission	RP09	12	120776.36	120800.36	120800.37	120896.83
01)Admission to Readmission	RP07	10	120781.05	120801.05	120801.06	120881.44
01)Admission to Readmission	RP05	8	120786.03	120802.03	120802.04	120866.35
01)Admission to Readmission	RP03	6	120800.94	120812.94	120812.95	120861.18
01)Admission to Readmission	RP02	5	120811.14	120821.14	120821.15	120861.34
01)Admission to Readmission	Generalized gamma	4	120965.59	120973.59	120973.59	121005.74
01)Admission to Readmission	Lognormal	3	121321.19	121327.19	121327.19	121351.31
01)Admission to Readmission	Generalized gamma (original)	4	121472.61	121480.61	121480.61	121512.76
01)Admission to Readmission	Log-logistic	3	121840.98	121846.98	121846.98	121871.10
01)Admission to Readmission	Weibull (AFT)	3	122215.75	122221.75	122221.76	122245.87
01)Admission to Readmission	Weibull (PH)	3	122215.75	122221.75	122221.76	122245.87
01)Admission to Readmission	Gamma	3	122337.40	122343.40	122343.40	122367.52
01)Admission to Readmission	Exponential	2	122716.58	122720.58	122720.59	122736.66
02)Readmission to Readmission2	Generalized F	5	36551.26	36561.26	36561.27	36595.07
02)Readmission to Readmission2	Generalized F(original)	5	36551.26	36561.26	36561.27	36595.07
02)Readmission to Readmission2	RP03	6	36550.53	36562.53	36562.55	36603.12
02)Readmission to Readmission2	RP04	7	36548.81	36562.81	36562.83	36610.15
02)Readmission to Readmission2	RP05	8	36547.29	36563.29	36563.31	36617.40
02)Readmission to Readmission2	RP08	11	36542.68	36564.68	36564.72	36639.08
02)Readmission to Readmission2	RP09	12	36540.81	36564.81	36564.85	36645.97
02)Readmission to Readmission2	RP06	9	36547.06	36565.06	36565.09	36625.93
02)Readmission to Readmission2	RP07	10	36545.25	36565.25	36565.29	36632.89
02)Readmission to Readmission2	RP10	13	36540.91	36566.91	36566.97	36654.84
02)Readmission to Readmission2	RP02	5	36560.67	36570.67	36570.68	36604.48
02)Readmission to Readmission2	Generalized gamma	4	36650.99	36658.99	36658.99	36686.04
02)Readmission to Readmission2	Lognormal	3	36675.82	36681.82	36681.82	36702.11
02)Readmission to Readmission2	Gompertz	3	36697.35	36703.35	36703.35	36723.64
02)Readmission to Readmission2	Generalized gamma (original)	4	36707.55	36715.55	36715.55	36742.60
02)Readmission to Readmission2	Log-logistic	3	36785.86	36791.86	36791.86	36812.15
02)Readmission to Readmission2	Weibull (AFT)	3	36905.33	36911.33	36911.33	36931.62
02)Readmission to Readmission2	Weibull (PH)	3	36905.33	36911.33	36911.33	36931.62
02)Readmission to Readmission2	Gamma	3	36930.62	36936.62	36936.62	36956.91
02)Readmission to Readmission2	Exponential	2	36978.44	36982.44	36982.44	36995.97
03)Readmission2 to Readmission3	Generalized F(original)	5	11508.11	11518.11	11518.14	11546.14
03)Readmission2 to Readmission3	Generalized F	5	11508.11	11518.11	11518.14	11546.14
03)Readmission2 to Readmission3	RP02	5	11511.21	11521.21	11521.24	11549.24
03)Readmission2 to Readmission3	RP03	6	11511.36	11523.36	11523.40	11557.00
03)Readmission2 to Readmission3	RP04	7	11510.56	11524.56	11524.62	11563.81
03)Readmission2 to Readmission3	RP05	8	11510.45	11526.45	11526.52	11571.30
03)Readmission2 to Readmission3	RP06	9	11509.57	11527.57	11527.66	11578.03
03)Readmission2 to Readmission3	RP07	10	11509.10	11529.10	11529.21	11585.17
03)Readmission2 to Readmission3	RP08	11	11507.51	11529.51	11529.64	11591.18
03)Readmission2 to Readmission3	RP09	12	11506.53	11530.53	11530.68	11597.80
03)Readmission2 to Readmission3	RP10	13	11505.29	11531.29	11531.47	11604.18
03)Readmission2 to Readmission3	Generalized gamma	4	11530.60	11538.60	11538.62	11561.02
03)Readmission2 to Readmission3	Lognormal	3	11540.98	11546.98	11546.99	11563.80
03)Readmission2 to Readmission3	Generalized gamma (original)	4	11552.49	11560.49	11560.51	11582.92
03)Readmission2 to Readmission3	Gompertz	3	11568.15	11574.15	11574.16	11590.97
03)Readmission2 to Readmission3	Log-logistic	3	11578.40	11584.40	11584.42	11601.22
03)Readmission2 to Readmission3	Weibull (AFT)	3	11621.90	11627.90	11627.91	11644.72
03)Readmission2 to Readmission3	Weibull (PH)	3	11621.90	11627.90	11627.91	11644.72
03)Readmission2 to Readmission3	Gamma	3	11626.88	11632.88	11632.90	11649.70
03)Readmission2 to Readmission3	Exponential	2	11630.63	11634.63	11634.63	11645.84
04)Readmission3 to Readmission4	Generalized F(original)	5	3607.24	3617.24	3617.34	3639.61
04)Readmission3 to Readmission4	RP02	5	3607.36	3617.36	3617.45	3639.72
04)Readmission3 to Readmission4	Generalized F	5	3607.61	3617.61	3617.70	3639.97
04)Readmission3 to Readmission4	RP03	6	3607.49	3619.49	3619.62	3646.32
04)Readmission3 to Readmission4	RP04	7	3607.08	3621.08	3621.26	3652.39
04)Readmission3 to Readmission4	RP05	8	3606.09	3622.09	3622.32	3657.87
04)Readmission3 to Readmission4	RP06	9	3605.59	3623.59	3623.88	3663.84
04)Readmission3 to Readmission4	RP07	10	3605.51	3625.51	3625.86	3670.24
04)Readmission3 to Readmission4	Lognormal	3	3619.70	3625.70	3625.74	3639.12
04)Readmission3 to Readmission4	RP08	11	3604.67	3626.67	3627.08	3675.86
04)Readmission3 to Readmission4	Generalized gamma	4	3619.30	3627.30	3627.37	3645.19
04)Readmission3 to Readmission4	RP09	12	3604.95	3628.95	3629.44	3682.61
04)Readmission3 to Readmission4	Generalized gamma (original)	4	3621.96	3629.96	3630.02	3647.85
04)Readmission3 to Readmission4	RP10	13	3605.14	3631.14	3631.71	3689.28
04)Readmission3 to Readmission4	Gompertz	3	3626.04	3632.04	3632.08	3645.46
04)Readmission3 to Readmission4	Log-logistic	3	3626.47	3632.47	3632.51	3645.89
04)Readmission3 to Readmission4	Exponential	2	3640.60	3644.60	3644.62	3653.54
04)Readmission3 to Readmission4	Weibull (AFT)	3	3639.23	3645.23	3645.27	3658.65
04)Readmission3 to Readmission4	Weibull (PH)	3	3639.23	3645.23	3645.27	3658.65
04)Readmission3 to Readmission4	Gamma	3	3640.18	3646.18	3646.22	3659.60

Session Info

Sys.getenv("R_LIBS_USER")

## [1] "C:/Users/CISS Fondecyt/OneDrive/Documentos/R/win-library/4.0"

rstudioapi::getSourceEditorContext()

## Document Context: 
## - id:        'EA74E2C1'
## - path:      'G:/Mi unidad/Alvacast/SISTRAT 2019 (github)/SUD_CL/Matching_Process15_JUN_21.Rmd'
## - contents:  <1528 rows>
## Document Selection:
## - [1476, 4] -- [1476, 4]: ''

#save.image("G:/Mi unidad/Alvacast/SISTRAT 2019 (github)/mult_state_carla.RData")

if (grepl("CISS Fondecyt",rstudioapi::getSourceEditorContext()$path)==T){
    save.image("C:/Users/CISS Fondecyt/OneDrive/Escritorio/SUD_CL/mult_state_15_jun.RData")
  } else if (grepl("andre",rstudioapi::getSourceEditorContext()$path)==T){
    save.image("C:/Users/andre/Desktop/SUD_CL/mult_state_15_jun.RData")
  } else if (grepl("E:",rstudioapi::getSourceEditorContext()$path)==T){
    save.image("E:/Mi unidad/Alvacast/SISTRAT 2019 (github)/mult_state_15_jun.RData")
  } else {
    save.image("G:/Mi unidad/Alvacast/SISTRAT 2019 (github)/mult_state_15_jun.RData")
  }

sessionInfo()

## R version 4.0.2 (2020-06-22)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19042)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=Spanish_Chile.1252  LC_CTYPE=Spanish_Chile.1252   
## [3] LC_MONETARY=Spanish_Chile.1252 LC_NUMERIC=C                  
## [5] LC_TIME=Spanish_Chile.1252    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] gurobi_9.1-0            slam_0.1-47             radiant.update_1.4.1   
##  [4] Rfast_2.0.3             RcppZiggurat_0.1.6      Metrics_0.1.4          
##  [7] muhaz_1.2.6.1           flexsurv_2.0            mstate_0.3.1           
## [10] DiagrammeR_1.0.6.1.9000 Amelia_1.7.6            Rcpp_1.0.5             
## [13] igraph_1.2.6            eha_2.8.1               cobalt_4.2.3           
## [16] MatchIt_3.0.2           tableone_0.12.0         stargazer_5.2.2        
## [19] reshape2_1.4.4          gridExtra_2.3           foreign_0.8-80         
## [22] survMisc_0.5.5          ggfortify_0.4.10        survminer_0.4.8        
## [25] ggpubr_0.4.0            epiR_1.0-15             forcats_0.5.0          
## [28] purrr_0.3.4             readr_1.3.1             tibble_3.0.3           
## [31] tidyverse_1.3.0         dplyr_1.0.1             treemapify_2.5.3       
## [34] sf_0.9-3                ggiraph_0.7.0           finalfit_1.0.1         
## [37] lsmeans_2.30-0          emmeans_1.4.8           RColorBrewer_1.1-2     
## [40] panelr_0.7.3            lme4_1.1-23             Matrix_1.2-18          
## [43] data.table_1.13.0       codebook_0.9.2          devtools_2.3.0         
## [46] usethis_1.6.1           sqldf_0.4-11            RSQLite_2.2.0          
## [49] gsubfn_0.7              proto_1.0.0             broom_0.7.0            
## [52] zoo_1.8-8               rbokeh_0.5.1            janitor_2.0.1          
## [55] plotly_4.9.2.1          kableExtra_1.1.0        Hmisc_4.4-0            
## [58] Formula_1.2-3           survival_3.1-12         lattice_0.20-41        
## [61] ggplot2_3.3.2           stringr_1.4.0           stringi_1.4.6          
## [64] tidyr_1.1.1             knitr_1.29              matrixStats_0.56.0     
## [67] boot_1.3-25            
## 
## loaded via a namespace (and not attached):
##   [1] estimability_1.3    coda_0.19-3         acepack_1.4.1      
##   [4] bit64_0.9-7         multcomp_1.4-13     rpart_4.1-15       
##   [7] generics_0.0.2      callr_3.4.3         TH.data_1.0-10     
##  [10] mice_3.11.0         ggfittext_0.9.0     chron_2.3-55       
##  [13] bit_1.1-15.2        webshot_0.5.2       xml2_1.3.2         
##  [16] lubridate_1.7.9     assertthat_0.2.1    xfun_0.16          
##  [19] hms_0.5.3           evaluate_0.14       fansi_0.4.1        
##  [22] dbplyr_1.4.4        readxl_1.3.1        km.ci_0.5-2        
##  [25] DBI_1.1.0           htmlwidgets_1.5.1   ellipsis_0.3.1     
##  [28] crosstalk_1.1.0.1   backports_1.1.7     insight_0.9.0      
##  [31] survey_4.0          vctrs_0.3.2         remotes_2.2.0      
##  [34] sjlabelled_1.1.6    abind_1.4-5         withr_2.2.0        
##  [37] pryr_0.1.4          checkmate_2.0.0     prettyunits_1.1.1  
##  [40] cluster_2.1.0       lazyeval_0.2.2      crayon_1.3.4       
##  [43] labeling_0.3        pkgconfig_2.0.3     units_0.6-6        
##  [46] nlme_3.1-148        pkgload_1.1.0       nnet_7.3-14        
##  [49] rlang_0.4.7         lifecycle_0.2.0     sandwich_2.5-1     
##  [52] modelr_0.1.8        cellranger_1.1.0    tcltk_4.0.2        
##  [55] rprojroot_1.3-2     KMsurv_0.1-5        carData_3.0-4      
##  [58] reprex_0.3.0        base64enc_0.1-3     processx_3.4.3     
##  [61] png_0.1-7           viridisLite_0.3.0   parameters_0.8.2   
##  [64] KernSmooth_2.23-17  visNetwork_2.0.9    pander_0.6.3       
##  [67] blob_1.2.1          classInt_0.4-3      jpeg_0.1-8.1       
##  [70] rstatix_0.6.0       ggeffects_0.15.1    ggsignif_0.6.0     
##  [73] scales_1.1.1        memoise_1.1.0       magrittr_1.5       
##  [76] plyr_1.8.6          hexbin_1.28.1       compiler_4.0.2     
##  [79] snakecase_0.11.0    cli_2.0.2           ps_1.3.3           
##  [82] htmlTable_2.0.1     MASS_7.3-51.6       tidyselect_1.1.0   
##  [85] highr_0.8           mitools_2.4         jtools_2.0.5       
##  [88] yaml_2.2.1          latticeExtra_0.6-29 grid_4.0.2         
##  [91] tools_4.0.2         parallel_4.0.2      rio_0.5.16         
##  [94] rstudioapi_0.11     uuid_0.1-4          gistr_0.5.0        
##  [97] farver_2.0.3        sjPlot_2.8.4        digest_0.6.25      
## [100] quadprog_1.5-8      car_3.0-8           performance_0.4.8  
## [103] httr_1.4.2          gdtools_0.2.2       effectsize_0.3.2   
## [106] sjstats_0.18.0      colorspace_1.4-1    rvest_0.3.6        
## [109] fs_1.5.0            splines_4.0.2       statmod_1.4.34     
## [112] sessioninfo_1.1.1   systemfonts_0.2.3   xtable_1.8-4       
## [115] jsonlite_1.7.0      nloptr_1.2.2.2      testthat_2.3.2     
## [118] R6_2.4.1            pillar_1.4.6        htmltools_0.5.0    
## [121] glue_1.4.1          minqa_1.2.4         deSolve_1.28       
## [124] class_7.3-17        codetools_0.2-16    maps_3.3.0         
## [127] pkgbuild_1.1.0      mvtnorm_1.1-1       numDeriv_2016.8-1.1
## [130] curl_4.3            BiasedUrn_1.07      zip_2.1.1          
## [133] openxlsx_4.1.5      rmarkdown_2.6       desc_1.2.0         
## [136] munsell_0.5.0       e1071_1.7-3         labelled_2.5.0     
## [139] sjmisc_2.8.5        haven_2.3.1         gtable_0.3.0       
## [142] bayestestR_0.7.2